function Lbeam
clc
%close all
addpath('../');
addpath('../../');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% read mesh
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
projDir = '/home/hbui/kratos_janosch';
%coarse (structure) mesh
projName = 'Lbeam.gid';
fn = [projDir '/' projName '/' 'Lbeam.msh'];
% %fine mesh at crack zone
% projName = 'Lbeam3.gid';
% fn = [projDir '/' projName '/' 'mesh3.msh'];

mg = fem2.gid.GidMeshGenerator(fn);
mesh = mg.getMesh2d();
mesh.setNumberOfFields(2);
for i=1:mesh.countNodes()
    mesh.getNode(i-1).setId(i);
end
%plotMesh(mesh);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% conditions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
numTimeStep = 7;
% 7 for Lbeam.msh
% 5 for mesh3.msh
ddu = -0.175;
du = numTimeStep*ddu;
meshUtil = fem2.MeshUtilities(mesh);
c1 = fem2.Constraint([0 0]);
c2 = fem2.Constraint([1 0]);
c2.setValue(1,du);
mesh.setConstraint(c1,meshUtil.seekNodesOnSurface(0,1,0));
meshUtil.seekNode([1 -1]).setConstraint(c2);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% element declaration
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
M = fem2.Model(mesh);
nf = mesh.countFaces();
E = 2;
nu = 0.3;
H = -1;
ft = 0.91;
C = C_plane_stress(E,nu);
for i=1:nf
    e = fem2.MatlabStructuralElement(mesh.getFace(i-1));
    M.addElement(e);
    %create struct to contain element state variables
    es = createElementContainer(C,H,ft);
    elements(i) = es;
end
%deactivate and activate elements
for i=1:nf
    c = computeFaceCenter(mesh.getFace(i-1));
    if (c(1)>=0 && c(1)<=0.6 && c(2)>=-0.6 && c(2)<=-0.4)
        elements(i).isActivated = 1;
    else
        elements(i).isActivated = 0;
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% time steps iteration
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
M.preAnalysis();
nDofs = mesh.countDofs();
r = zeros(nDofs,1);
M.assembleR(r);
u = zeros(nDofs,1);
tol = 1e-6;
lambda = 0;
dlambda = 1/numTimeStep;
for i=1:numTimeStep
    lambda = lambda+dlambda;
    lr = lambda*r;
    %update model's states
    M.setX(lambda,u);
    elements = InitializeSolutionStep(elements,M);
    
    %assemble ri
    ri = assembleRi(elements,M);
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Newton-Raphson loop
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%
    cnt = 0;
    while 1
        cnt = cnt+1;
        disp(['substep ' num2str(cnt) '--------------']);
        %assemble kt
        kt = assembleKt(elements,M);
        
        du = kt\(lr-ri);
        u = u + du;
        M.setX(lambda,u);
        
        %update model's states
        elements = CalculateMaterialResponse(elements,M);
        
        %assemble ri
        ri = assembleRi(elements,M);
        
        %compute convergent criteria
        %res = norm(du)/norm(u);
        res = norm(r-ri);
        
        disp(['substep ' num2str(cnt) ' complete']);
        if res < tol
            disp(['Newton Raphson iteration converged in ' ...
                num2str(cnt) ' loops']);
            break;
        end
        num = 50;
        if cnt>num
            disp(['Newton Raphson iteration did not converge in ' ...
                num2str(num) ' loops']);
            break;
        end
    end
    
    elements = FinalizeSolutionStep(elements,M);

    findMaximumPrincipalStress(elements);
    
    responsePlot.x(i) = i*ddu;
    responsePlot.y(i) = elements(398).stress(1);
    
    disp(['time step ' num2str(i) ...
        ' finished with lambda = ' num2str(lambda)]);
end

figure
plot(responsePlot.x,responsePlot.y);

figure
plotMeshWithCrack(M,elements);

rmpath('../');
rmpath('../../');

